#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gmp.h>
#include <time.h>
#include <pthread.h>
#define ITERACOES 25
#define DEZMI 10000000

typedef struct{
    mpf_t a;
    mpf_t b;
    mpf_t t;
    mpf_t p;
}Data;

pthread_t threadA, threadB, threadT, threadP;
Data D[ITERACOES];


void *calcula_A(void *k) {
	int n = (int) k;
	mpf_t temp;
	mpf_init_set_ui(temp, 0L); //temp = 0
	
	//Value for A(n+1)
	mpf_add(temp, D[n].a, D[n].b);
    mpf_div_ui(D[n+1].a, temp, 2L);
	pthread_exit((void *) k);
}

void *calcula_B(void *k){
	int n = (int) k;
	mpf_t temp;
	mpf_init_set_ui(temp, 0L); //temp = 0
	
	//Value for B(n+1)
    mpf_mul(temp, D[n].a, D[n].b);
    mpf_sqrt(D[n+1].b, temp);
	pthread_exit((void *) k);
}
	

void *calcula_T(void *k){
	int n = (int) k;		
	mpf_t temp;
	mpf_init_set_ui(temp, 0L); //temp = 0
 
	//value for T(n+1)
    mpf_sub(temp, D[n].a, D[n+1].a);
    mpf_pow_ui(temp, temp, 2L);
    mpf_mul(temp, D[n].p, temp);
    mpf_sub(D[n+1].t, D[n].t, temp);
	pthread_exit((void *) k);
}

void *calcula_P(void *k){
	int n = (int) k;
	mpf_t temp;
	mpf_init_set_ui(temp, 0L); //temp = 0
	
	//value for P(n+1)
    mpf_mul_ui(D[n+1].p, D[n].p, 2L);
	pthread_exit((void *) k);
}


int main(){
    time_t start, end;
    int i;
    int n;
	int flag;
				
    mpf_set_default_prec((DEZMI*4)+16);
    mpf_t PI;
	mpf_t temp;
	

 	
    for(i=0; i<ITERACOES; i++){
        //mpf_init(D[i].a);
        mpf_init_set_d(D[i].a, 0.0);

        //mpf_init(D[i].b);
        mpf_init_set_d(D[i].b, 0.0);

        //mpf_init(D[i].t);
        mpf_init_set_d(D[i].t, 0.0);

        //mpf_init(D[i].p);
        mpf_init_set_d(D[i].p, 0.0);
    }

    //mpf_init(PI);
    //mpf_init(temp);
    mpf_init_set_ui(PI, 0L); // pi =0
    mpf_init_set_ui(temp, 0L); //temp = 0

    //Seta valores iniciais (a0, b0, t0 e p0)
    mpf_init_set_ui(D[0].a, 1L);

    mpf_sqrt_ui(temp, 2L); // temp = sqrt(2)
    mpf_ui_div(D[0].b, 1L, temp); // b0 = 1/sqrt(2)

    mpf_init_set_d(D[0].t, 0.25); //t0 = 0.25

    mpf_init_set_ui(D[0].p, 1L); // p0 = 1

    /*And let's the magic begin: */
	
	//Setar valores da rodada 1 de A,B e P
	 
	//Value for A(1)
    mpf_add(temp, D[0].a, D[0].b);
    mpf_div_ui(D[1].a, temp, 2L);

    //Value for B(1)
    mpf_mul(temp, D[0].a, D[0].b);
    mpf_sqrt(D[1].b, temp);


    //value for P(1)
    mpf_mul_ui(D[1].p, D[0].p, 2L);
	

    for(n=1; n<ITERACOES-1; n++){
       
        flag = pthread_create(&threadA, NULL, calcula_A, (void *) n);

        if(flag){
            printf("ERROR\n");
            exit(1);
        }
        flag = pthread_create(&threadB, NULL, calcula_B, (void *) n);

        if(flag){
            printf("ERROR\n");
            exit(1);
        }
        flag = pthread_create(&threadT, NULL, calcula_T, (void *) (n - 1));

        if(flag){
            printf("ERROR\n");
            exit(1);
        }
	
        flag = pthread_create(&threadP, NULL, calcula_P, (void *) n);

        if(flag){
            printf("ERROR\n");
            exit(1);
        }

    pthread_join(threadA, NULL);  
	pthread_join(threadB, NULL);
	pthread_join(threadT, NULL);
	pthread_join(threadP, NULL); 

        
    }


	//value for T(n+1)
    mpf_sub(temp, D[n-1].a, D[n].a);
    mpf_pow_ui(temp, temp, 2L);
    mpf_mul(temp, D[n-1].p, temp);
    mpf_sub(D[n].t, D[n-1].t, temp);


    //value for PI
    mpf_add(temp, D[n].a, D[n].b);
    mpf_pow_ui(PI, temp, 2L);
    mpf_mul_ui(temp, D[n].t, 4L);
    mpf_div(PI, PI, temp);
	
	
	mpf_out_str( stdout, 10, (DEZMI*4)+16, PI );
    		
    

    return 0;
}
